%Name:          Airy_disc.m
%Last Modified: 23/02/2012
%Description:   Script M-file to compute the Fraunhofer diffraction pattern
%               from a circular aperture (Airy disc)
%Written by:    Mike Shaw
%Notes:         

function [x,PSF] = Airy_disc(p,alpha)

lambda = 488e-9;    %wavelength of light, m
k = 2*pi/lambda;    %wavenumber of light, 1/m

x = -100*p:p:100*p;
[X,Y] = meshgrid(x,x);
[~,r] = cart2pol(X,Y);

N = numel(r);

l = k*tan(alpha)*abs(r);
PSF = zeros(numel(x),numel(x));     %vector of intensity values at image plane

for i=1:N
    if r(i)==0
    PSF(i) = 1;
    else
    PSF(i) = (2*besselj(1,l(i))/l(i))^2;
    end
end